## Loading required package: assertthat
## 
## Attaching package: 'assertthat'
## The following object is masked from 'package:tibble':
## 
##     has_name
## Loading required package: glue
## 
## Attaching package: 'glue'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## Loading required package: rlang
## Warning: package 'rlang' was built under R version 4.0.3
## 
## Attaching package: 'rlang'
## The following object is masked from 'package:assertthat':
## 
##     has_name
## The following objects are masked from 'package:purrr':
## 
##     %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
##     flatten_lgl, flatten_raw, invoke, list_along, modify, prepend,
##     splice
## 
## -- Column specification --------------------------------------------------------
## cols(
##   id = col_double(),
##   code = col_character(),
##   naam_nederlands = col_character(),
##   naam_wetenschappelijk = col_character(),
##   taxn_vrs_key = col_character(),
##   taxongroep = col_character(),
##   usageKey = col_double(),
##   scientificName = col_character(),
##   rank = col_character(),
##   order = col_character(),
##   matchType = col_character(),
##   phylum = col_character(),
##   kingdom = col_character(),
##   genus = col_character(),
##   class = col_character(),
##   confidence = col_double(),
##   synonym = col_logical(),
##   status = col_character(),
##   family = col_character()
## )
## i Using ',' as decimal and '.' as grouping mark. Use `read_delim()` for more control.
## 
## -- Column specification --------------------------------------------------------
## cols(
##   substraat = col_character(),
##   substrate = col_character()
## )

Exploratory data analysis

## `summarise()` ungrouping output (override with `.groups` argument)
phylum n_obs n_soorten n_hokken
Anthocerotophyta 218 4 123
Bryophyta 95349 410 796
Marchantiophyta 17462 109 727 

## `summarise()` regrouping output by 'phylum' (override with `.groups` argument)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Alternative Telfer method

Goed onderzocht = 40 of meer soorten in een uurhok per periode

Uurhok moet in beide perioden goed onderzocht zijn.

## Joining, by = "substraat"
phylum substrate number of species
Bryophyta artificial stone 84
Bryophyta bryophyte 3
Bryophyta decaying vegetation 21
Bryophyta decorticated wood 48
Bryophyta epiphytic on living wood 73
Bryophyta floating on water 1
Bryophyta gravel or sand 106
Bryophyta mineral soil 210
Bryophyta peat 46
Bryophyta rock, hard 95
Bryophyta rock, soft 38
Bryophyta soil on rock 68
Marchantiophyta artificial stone 9
Marchantiophyta bryophyte 17
Marchantiophyta decaying vegetation 18
Marchantiophyta decorticated wood 16
Marchantiophyta epiphytic on living wood 15
Marchantiophyta floating on water 3
Marchantiophyta gravel or sand 34
Marchantiophyta mineral soil 67
Marchantiophyta peat 30
Marchantiophyta rock, hard 14
Marchantiophyta rock, soft 14
Marchantiophyta soil on rock 20
phylum ell_l number of species
Bryophyta 1 9
Bryophyta 2 19
Bryophyta 3 37
Bryophyta 4 86
Bryophyta 5 125
Bryophyta 6 188
Bryophyta 7 188
Bryophyta 8 126
Bryophyta 9 15
Marchantiophyta 1 1
Marchantiophyta 3 20
Marchantiophyta 4 51
Marchantiophyta 5 51
Marchantiophyta 6 60
Marchantiophyta 7 56
Marchantiophyta 8 18
phylum ell_r number of species
Bryophyta 1 10
Bryophyta 2 63
Bryophyta 3 88
Bryophyta 4 67
Bryophyta 5 95
Bryophyta 6 138
Bryophyta 7 214
Bryophyta 8 101
Bryophyta 9 17
Marchantiophyta 1 28
Marchantiophyta 2 60
Marchantiophyta 3 20
Marchantiophyta 4 52
Marchantiophyta 5 30
Marchantiophyta 6 38
Marchantiophyta 7 21
Marchantiophyta 8 6
Marchantiophyta 9 2
phylum ell_f number of species
Bryophyta 1 12
Bryophyta 2 12
Bryophyta 3 61
Bryophyta 4 152
Bryophyta 5 187
Bryophyta 6 176
Bryophyta 7 50
Bryophyta 8 69
Bryophyta 9 51
Bryophyta 10 14
Bryophyta 11 4
Bryophyta 12 5
Marchantiophyta 4 15
Marchantiophyta 5 51
Marchantiophyta 6 52
Marchantiophyta 7 50
Marchantiophyta 8 47
Marchantiophyta 9 32
Marchantiophyta 10 9
Marchantiophyta 11 1
phylum ell_n number of species
Bryophyta 1 32
Bryophyta 2 192
Bryophyta 3 139
Bryophyta 4 195
Bryophyta 5 137
Bryophyta 6 71
Bryophyta 7 27
Marchantiophyta 1 43
Marchantiophyta 2 97
Marchantiophyta 3 45
Marchantiophyta 4 36
Marchantiophyta 5 23
Marchantiophyta 6 8
Marchantiophyta 7 5
phylum ell_t number of species
Bryophyta 2 48
Bryophyta 3 269
Bryophyta 4 142
Bryophyta 5 155
Bryophyta 6 108
Bryophyta 7 58
Bryophyta 8 13
Marchantiophyta 2 27
Marchantiophyta 3 113
Marchantiophyta 4 51
Marchantiophyta 5 43
Marchantiophyta 6 13
Marchantiophyta 7 2
Marchantiophyta 8 8

Only remove substrate levels with less than 5 species (for Ellenberg variables this is not necessary because they are included as continuous covariate).

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

phylum substrate n
Bryophyta artificial stone 84
Bryophyta decaying vegetation 21
Bryophyta decorticated wood 48
Bryophyta epiphytic on living wood 73
Bryophyta gravel or sand 106
Bryophyta mineral soil 210
Bryophyta peat 46
Bryophyta rock, hard 95
Bryophyta rock, soft 38
Bryophyta soil on rock 68
Marchantiophyta artificial stone 9
Marchantiophyta bryophyte 17
Marchantiophyta decaying vegetation 18
Marchantiophyta decorticated wood 16
Marchantiophyta epiphytic on living wood 15
Marchantiophyta gravel or sand 34
Marchantiophyta mineral soil 67
Marchantiophyta peat 30
Marchantiophyta rock, hard 14
Marchantiophyta rock, soft 14
Marchantiophyta soil on rock 20 
## Warning: Ignoring unknown aesthetics: text
## `geom_smooth()` using formula 'y ~ x'
## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: Ignoring unknown aesthetics: text
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 7 rows containing non-finite values (stat_smooth).
## # A tibble: 124 x 2
##    n_hok_1980_1999     n
##              <int> <int>
##  1               0   241
##  2               5    35
##  3               6    17
##  4               7    18
##  5               8    12
##  6               9    22
##  7              10    14
##  8              11    12
##  9              12    15
## 10              13     4
## # ... with 114 more rows
## # A tibble: 145 x 2
##    n_hok_2000_2019     n
##              <int> <int>
##  1               0    18
##  2               1    86
##  3               2    63
##  4               3    33
##  5               4    39
##  6               5    20
##  7               6    22
##  8               7    14
##  9               8    25
## 10               9    16
## # ... with 135 more rows

Keeping only species that occur in both periods (much improves model validation). (Alternative could be to include a factor with levels for lost, persisting, gained in the model to explain lost and gained species, but maybe see this as a separate model with the full dataset!)

Overdispersion (residual deviance divided by residual degrees of freedom >> 1)! Quasi-binomial or beta-binomial model needed.

Should we include interaction between baseline proportion and substrate or ellenberg variables?

Model formulation:

\[ g(\mu_{i,s2}) = (\beta_0 + \beta_1\log\frac{\mu_{i,s1}}{1-\mu_{i,s1}} + \beta_L\text{ell}_L + \dots)*\text{phylum}_k + \mathbf{b}_{o,j} + \mathbf{b}_{1,j}\text{phylum}_k \] where:

\[ \begin{aligned} \mu_{i,s2} \equiv \mathbf{E}(Y_{i,s2}|\text{trials},\mathbf{b})&& \text{expected value}\\ g(\mu_{i,s2}) = \log\frac{\mu_{i,s2}}{1-\mu_{i,s2}} && \text{logit-link} \\ Y_{i,s2} \sim \text{BetaBinomial}(\mu_{i,s2}, \phi_{i,s2}) && \text{Beta-Binomial distribution}\\ \begin{bmatrix}\mathbf{b}_{o,j}\\ \mathbf{b}_{1,j}\end{bmatrix}\sim N(0,\Sigma) && \text{random intercept and slope}\\ \Sigma = \begin{bmatrix}\sigma^2_{0,j}&\sigma_{01,j}\\\sigma_{01,j}&\sigma^2_{1,j}\end{bmatrix} && \text{variance-covariance matrix} \end{aligned} \]

After rearranging (for one phylum) and exponentiating the linear predictor, ignoring all random effects, we obtain the following odds-ratio:

\[ \frac{Odds_{\mu_{i,s2}}}{Odds^{\beta_1}_{\mu_{i,s1}}} = \exp(\beta_0 + \beta_L\text{ell}_L + \dots) \] which allows us to predict the effect of Ellenberg values on the

## Registered S3 method overwritten by 'broom.mixed':
##   method      from 
##   tidy.gamlss broom
## Warning in tidy.brmsfit(alt_telfer): some parameter names contain underscores:
## term naming may be unreliable!
effect component group term estimate std.error conf.low conf.high
fixed cond NA (Intercept) 0.6067183 0.1807537 0.2571887 0.9683907
fixed cond NA qlogisn_hok_1980_1999Dn_hok_tot_1980_1999 0.8910590 0.0172844 0.8564729 0.9235065
fixed cond NA phylumMarchantiophyta -0.7128882 0.3364575 -1.3682412 -0.0538964
fixed cond NA ell_l -0.0205340 0.0173167 -0.0537886 0.0123278
fixed cond NA ell_n -0.0200793 0.0199998 -0.0590172 0.0188560
fixed cond NA ell_f -0.0559915 0.0162536 -0.0883449 -0.0248281
fixed cond NA ell_t 0.0320731 0.0212892 -0.0099047 0.0728306
fixed cond NA qlogisn_hok_1980_1999Dn_hok_tot_1980_1999:phylumMarchantiophyta 0.0146853 0.0427090 -0.0694389 0.0973873
fixed cond NA phylumMarchantiophyta:ell_l -0.0131476 0.0359741 -0.0843111 0.0553381
fixed cond NA phylumMarchantiophyta:ell_n 0.0505936 0.0435477 -0.0349141 0.1362271
fixed cond NA phylumMarchantiophyta:ell_f -0.0023349 0.0375804 -0.0759480 0.0726301
fixed cond NA phylumMarchantiophyta:ell_t 0.0609244 0.0464215 -0.0318811 0.1528604
ran_pars cond substrate sd__(Intercept) 0.1596420 0.0535199 0.0790869 0.2866274
ran_pars cond substrate sd__phylumMarchantiophyta 0.1013226 0.0809064 0.0035919 0.2966237
ran_pars cond substrate cor__(Intercept).phylumMarchantiophyta 0.2972960 0.5161680 -0.8468412 0.9757969

From Telfer et al (2002) \(\text{logit}(P_2k) = a + b \text{logit}(P_1k) + Z_k\): > The intercept a measures change in recording probability(on a logit scale) and the slope b measures the extent to which that change depends on the recording probability in the first survey. Parameters a and b may include the effects both of biological change and of change in recorder behaviour.

## Adding missing grouping variables: `phylum`
## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(ellenberg)` instead of `ellenberg` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## Adding missing grouping variables: `phylum`
## Adding missing grouping variables: `phylum`
## Adding missing grouping variables: `phylum`